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ABSTRACT 

The bulk of the solar chromosphere is weakly ionized and interactions between 
ionized particles and neutral particles likely have significant consequences for the 
thermodynamics of the chromospheric plasma. We investigate the importance 
of introducing neutral particles into the MHD equations using numerical 2.5D 
radiative MHD simulations obtained with the Bifrost code. The models span 
the solar atmosphere from the upper layers of the convection zone to the low 
corona, and solve the full MHD equations with non-grey and non-LTE radiative 
transfer, and thermal conduction along the magnetic field. The effects of partial 
ionization are implemented using the generalized Ohm's law, i.e., we consider the 
eflFects of the Hall term and ambipolar diflFusion in the induction equation. The 
approximations required in going from three fiuids to the generalized Ohm's law 
are tested in our simulations. The Ohmic diffusion, the Hall term, and ambipolar 
diffusion show strong variations in the chromosphere. These strong variations of 
the various magnetic diffusivities are absent or significantly underestimated when, 
as has been common for these types of studies, using the semi-empirical VAL-C 
model as a basis for estimates. In addition, we find that differences in estimating 
the magnitude of ambipolar diffusion arise depending on which method is used to 
calculate the ion-neutral collision frequency. These differences cause uncertainties 
in the different magnetic diffusivity terms. In the chromosphere, we find that the 
ambipolar diffusion is of the same order of magnitude or even larger than the 
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numerical diffusion used to stabilize our code. As a consequence, ambipolar 
diffusion produces a strong impact on the modeled atmosphere. Perhaps more 
importantly, it suggests that at least in the chromospheric domain, self-consistent 
simulations of the solar atmosphere driven by magneto-convection can accurately 
describe the impact of the dominant form of resistivity, i.e., ambipolar diffusion. 
This suggests that such simulations may be more realistic in their approach to the 
lower solar atmosphere (which directly drives the coronal volume) than previously 
assumed. 

Subject headings: Magnetohydrodynamics MHD — Methods: numerical — Ra- 
diative transfer — Sun: atmosphere — Sun: magnetic field 



Introduction 



Most of the models and simulations of the solar atmosphere solve the magnetohy- 
drodynamics (MHD) equations, implicitly assuming that the plasma is magnetized, i.e., 
fully ionized or with the ion-neutral collision frequency lower than the ion gyrofrequency 



( ISchaffenberger et al 



2005 



Vogler et al.ll2005l : IStein fc Nordlundll2006l : iGudiksen et al 



2011 



among others). However, since the photosphere and parts of the chromosphere are unmag- 
netized, i.e., the ions are not necessarily tied to the field lines, we expect that neutral 



particles can have a s ignific ant impact on the dynamics of this region (j Vernazza et al.lll981 



Fontenla et al.l Il990l . 1 19931 ). Therefore, it is likely that, under some conditions, the pho- 
tosphere and chromosphere should be treated as a three component fluid, where the dy- 
namics of the neutrals, ions, and electrons are treated separately. Under the assumption 
of a weakly ionized plasma, one can return to a one-component fluid. However, ne w terms 
known as the Hall term and ambipo lar diffusion appear in the induction equation ([Parker 
19631 . 120071 : IPandey fc Wardld 120081 ) . The latter term is a cons equence of the ion-neutral 



dissipation which can be derived from the Cowling resistivity ( iKhodachenko et al. 



2006; 


Leake & Arber 


2006) 


Ohm's law ( 


CowlingI 


1957) 



2004 



A large number of papers in r ecent years have inves tigated the effects of the ion- neutral 
interactions on single fluid MHD. iLeake fc Arberl ( 120061 ) simulated 2.5D simulations of flux 
emergence and observed that the ambipolar diffusion leads to an increase of the rates of 
magnetic field emergence and a resultant magnetic field that is much more diffuse than 
the case with only Ohmic diffusivity. In addition, the magnetic field that emerges into the 
corona is found to be more force-free, since currents are aligned to the field. Th is is because 



ambipolar diffusion acts on the currents perpendicular to the magnetic field. lArber et al. 
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( 120071 ) extended this simulation to 3D where the previous results were confirmed, and in 
addition found that, as a result of including neutrals, fiux emergence lifts less chromospheric 
material to great heights. This effect suppresses the Ray leigh- Taylor instability between the 
emerging fiux and the corona. 

The interaction between ions and neutrals can also dissipate Alfven waves as result of the 
small but finite coupling time between ions and neutrals. This type of damping can heat and 
accel e rate the plasma in the upper chromosphe r e and in spicules (IDe Pontieu fc Haerendel 



19981 : iDe Pontieul Il999l : IJames fc Erdelvil l2002l : Ijames et all 12003 



Erdelvi fc James 



2om 



and incur wave ene rgy leakage at t h e footpoints of coronal loo ps ( iDe Pontieu et al.ll200ll ) 
and in the network (jGoodmanll2000l ) . iKhodachenko et al.l (120041 ). using the temperature and 
density structure from the ID VAL-C model, concluded that the collisional friction damping 
of MHD waves is often more important than the viscous damping for waves propagating 
in the partially ionized plasmas of the solar photosphere, chro mosphere and prominences. 
Estimates of the efiiciency of the damping of waves were made by iLeake et al.l ( 120051 ) as well 
as by the previous authors. 



Pandey fc Wardld (120081 ) determine that waves can be affected by the Hall term at 



both low and high fractional ionization, because the Hall regime wave damping is inversely 
proportional to the fractional ionization. Thus Hall term may also be important at high frac- 
tional ionization in contrast to ambipolar diffusion which is important only at low fractional 
ionization. 



Khomenko fc CoUadosI ( 120121 ) performed various simplified scenarios where they studied 



the impact of the ambipolar diffusion in the chromosphere. They conclude that current 
dissipation enhanced by the action of ambipolar diffusion is an important process that is 
able to provide a significant energy input into the chromosphere. Heating from ambipolar 
diffusion leads to thermodynamic evolution in the chromosphere on timescales of about 10- 
100 seconds. 

All the models above, even the 2D and 3D models, are based on a ID semi-empirical 
atmosphere (e.g., VAL-C), and/or a simplified approach to the energy balance in the chromo- 
sphere (adiabatic or Newtonian cooling). In addition, none of the partial ionization effects 
have been considered in full magneto-convection simulations. ICheung fc CameronI (120121 ) 
have made progress in this direction and performed full magneto-convection simulations of 
an umbra taking into account partial-ionization effects. However, their simulations only 
extend up to the upper photosphere. 



In this paper, we use the Bifrost code ( IGudiksen et al.ll201ll ) to create a self-consistent 
and fully dynamic model atmosphere of the sun, from the convection zone to the corona, to 
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consider the importance of the HaU term and ambipolar diffusion relative to the Ohmic and 
artificial diflFusion. Unlike other models, Bifrost includes an advanced treatment of radiative 
losses in the chromosphere based on recipes derived from dynamic non-LTE radiative ID 
hydrodynamic simulations. Such a treatment is crucial for a consideration of the effects 
of partial ionization, as shown in what follows. The code and the implementation of the 
generalized Ohm's law are described in Section [21 The tests performed for the code valida- 
tion are discussed in Section 12. 2[ We describe the different forms of diffusion in 2D MHD 
simulations in Section 13. 1[ Fin ally, the various simplifica tions made in order to obtain the 
generalized Ohm's law following IPandey fc Wardld ( 120081 ) have been investigated and tested 
for the 2D MHD simulations in Section 13. 2[ The paper finishes by addressing the conclusions 
and discussion. 



2. Equations and numerical method 



The magnetic upper-photosphere and chromosphere is weakly ionized and the interac- 
tion between ionized part icles and neutral part icles potentially has important consequences 
for the thermodynamics (IFontenla et al.lll993l ) of this region. We investigate these conse- 
quences in the solar atmosphere. In order to model the solar atmosphere we solve the MHD 
equations in 2.5D. The model spans from the upper layers of the convection zone to the low 
corona. We have implemented the effects of partial ionization into the induction equation 
through the Hall and ambipolar diffusion terms as described below. 



The Bifrost ( IGudiksen et al.ll201ll ) code is a staggered mesh, explicit code that solves the 
MHD partial differential equations, including non-LTE and non-grey radiative transfer with 
scattering, and conduction along the magnetic field lines. A lookup table, based on LTE, 
is used to compute the temperature, pressure, opacities and other radiation quantities, and 
ionization state, given the pressure and the internal energy of the plasma. Spatial derivatives 
and the interpolation of variables are done using high order polynomials. The equations are 
ste pped forward i n tim e using the explicit third order predictor-corrector procedure described 
by iHyman et al.l ( 119791 ). In order to suppress numerical noise, high-order artificial diffusion 



is added both in the forms of a viscosity and in the form of a magnetic diffusivity (see 
Gudiksen et al.l 1201 iL for details). 



The Bifrost code includes an advanced treatment of the effects of radiation on the local 
energy balance, which is crucial if one wants to accurately determine the ionization degree. 
The radiative fiux divergence from the photosphere and lower chromosphere is obtained by 
angle and wavelength integration of the transport equation assuming isotropic opacities and 
emissivities. The transport equation assumes that opacities are in LTE using four group 
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mean opacities to cover the entire spectrum (INordlundl 1 19821 ) . This is done by formulating 
the transfer equation for each of the four bins, calculating a mean source function in each 
bin. These source functions contain an approximate coherent scattering term and an exact 
contribution from thermal emissivity. The resulting 3D scattering problems are solved by 
iteration, based on o ne-ray approxim ation in the angle integral for the mean intensity, a 
method developed bv ISkartlienI ( 120001 ). 



In the mid and upper chromosphere, the Bifrost code includes non-LTE radiative losses 
from tabulated hydrogen continua, hydrogen lines, and lines from singly ionized calcium as 
functions of temperature and column mass (jCarlsson fc LeenaartsI l2012l ). These radiative 
losses depend on the computed non-LTE escape probability as a function of column mass 
and are based on a ID dynamical chromospheric model in which the radiative losses are 
computed in detail dCarlsson fc Steinlll992[ Il994l . Il997l . l2002h . 



The energy dissipated by Joule heating is given by Qjouie = E • J where the electric field 
E is calculated from the current J, taking into account hi^h-order artificial re sistivity. The 
resistivity is computed using a hyper-diffusion operator ( IGudiksen et al.ll201ll ). This entails 
that the Joule heating due artificial diffusion is set proportional to the current squared times 
a factor that becomes large (of order 10) when magnetic field gradients are large, and is 
unity otherwise. 



2.1. Generalized Ohm's Law theory 

2 J J. Multi-fluid 

Most codes treat the solar atmosphere as a single fiuid where collisional frequencies 
are considered sufficient to ensure that all species are well coupled and that the momentum 
and energy equations can be added without the introduction of frictional terms or similar. 
However, as chromosp heric temperatures are likely to drop to a few 10^ K or even lower 
(ILeenaarts et al.ll201ll ). there is a high probability that plasma is only partially ionized and 
that "slippage" effects could become important. In this case the MHD equations should 
be treated by considering the plasma to consist of three fiuids: ions, electrons and neutral 
particles. The mass density for each type of particle is governed by the continuity equation 
applied to each species separately: 



|^ + V.ftuj=0 



(1) 
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where pj = rrijUj^ Uj, rij and rrij are the mass density, velocity, number density and particle 
mass of the ion, electron, and neutral species, i.e., j = i, e, n respectively. The mass transfer 
term as result of ionization and recombination has been neglected. This approximation is 
valid for a one-fluid approach if the system is in ionization balance, and there is no decoupling 
of ions and neutrals. 

The momentum equation, written in SI units, for each species, is as follows: 



Pi + Ui • V j Ui = -\/Pi + riiZqe (E + Ui x B) - ^ z/^^(ui - Uj) 

^ ^ j=e,n 

Pe + Ue • V J Ue = "VPg " riiQe (E + Ue X B) - Pe Z^ej(Ue - Uj) 

^ ^ j=i,n 

Pn + Un • Un = "VPn - Pn ^nj{Un - Uj) 



(2) 
(3) 
(4) 



where and Z are respectively the electron charge and ion charge. E and B are the electric 
and magnetic fleld and Pj = UikTi is the partial pressure of the jth species, k is Boltzmann's 
constant, and Vij is the collision frequency for species i with species j. We assume that 
collisions are sufliciently numerous that the ion and electron temperature can be considered 
the same (T^ = Tg). All three equations are linked through the last term, i.e., the exchange 
of momentum between the particles, where we have ignored the thermal force. In a similar 
manner as for the continuity equation, the momentum transfer term as result of ionization 
and recombination has been neglected. 

The number of equations thus increases considerably c ompared with s in gle fluid MHD, 
but by considerin g som e simpliflcations, as described by ICowlingI ( 119571 ): IParkerl ( 120071 ): 
Pandey fc Wardld ( 120081 ). one can easily generalize the MHD equations for each species to 
a single fluid (see below too). Therefore, the mass density is governed by the continuity 
equation for the bulk fluid as follows: 



^ + V-pu = (5) 

where the density for the bulk fluid is the sum of the different particle densities {p = Pi + 
Pe + Pn), and considering pi/p » Pe/p^ then p ^ Pi + Pn- In a similar manner, the velocity 
of the bulk fluid is u {pi^i + Pn^n)/p^ where the electron inertia is implicitly neglected in 
the deflnition of the bulk velocity. If we deflne the neutral density fraction {D = Pn/p)^ then 
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u ^ (1 — D)ui + Dun. Finally, the current density is given by J = rieqei^i — Ue) (assuming 
singly charged ions). Since Equation [5] is the same as for the single fluid formulation, the 
continuity equation does not need any modiflcation in the Bifrost code. 



Following Pandey fc Wardld (120081 ) , the single fluid momentum equation can be recov- 
ered if we neglect the eflFects of the electron inertia. Because it is implicitly neglected in 
the deflnition of the bulk velocity, it can also be neglected in the continuity and momentum 
equations. For simplicity, the ions are assumed to be singly charged, and we adopt charge 
neutrality (n^ = ng). In addition, the drift momentum is assumed to be considerably smaller 
than the fast momentum {^P\/v\ + c^) so that: 

Pipnu]^ «p^(^^ + c^), (6) 



where Ud = Ui — Un, = ^1 yJPoP^ and = ^J'^Pj p are respectively the drift, Alfven and 
sound velocities in the bulk fluid; Po is the vacuum permeability, and 7 is the ratio of specific 
heats. When the plasma does not fulfill Equation [6l the fiuids are strongly decoupled. This 
happens when the ion-neutral collision frequency is low. When the drift momentum is low, 
the drift momentum can be neglected for small dynamical frequencies (i.e., changes of the 
plasma properties on timescales commensurate with such frequencies): 



^ - , . , ^ni (7) 



where /3e = the ratio of the cyclotron frequency and the coUisional frequency. With 
these assumptions w e recover the single fiuid momentum equation as it is implemented in 



the Bifrost code (see lPandey fc Wardld l2008l . for details): 



p ( ^ + u • V ) u = - VP + J X B (8) 



2.1.2. Induction equation 
The Ohmic diffusion. Hall term, and ambipolar diffusion are given by 



Tiohm = - (9) 
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\B\ 



Vhall — 



\B\Pnlpf mPn/p? 



Vamb 



Pn^n 



The electrical conductivity (a) in the absence of a magnetic field is 



where the sums of the collision frequencies are written 



(10) 

(11) 



(12) 



(13) 
(14) 



In order to obtain the induction equation the following assumptions are made: 



First, the electric field 



E + Ui X B 



VPe J J X B me^en 



+ - + 



-Ud 



(15) 



is deduced from the electron momentum equation assuming zero electron inertia and 
is expressed in the ion's rest frame. 



• The plasma obeys 



• The term 



peiyen « Pi^iw 



PiPn 



(ud • V)ui - (ui • V)ud 



can be neglected when the dynamical frequency of the plasma is small 

^ < l^niP/Pi 



(16) 



(17) 



(18) 



-9- 



• Biermann's battery contribution, from the VPg/^e^e term in Equation [T5l is neglected. 

• The term D(5i/ (3^ is smaU and of order < 10~^, so it is also neglected, where (5j = uJcj/iyj 
is the ratio of the cyclotron frequency to the sum of the jth particle collision frequency. 

• Finally, terms due to the pressure gradient VP x B are negligible compared to the 
induction term u x B when the dynamical frequency is small: 

^ < (^] (19) 
Under these assumptions the electric field is defined as 



_ J J X B 

E = - + 



,J X B X B 



(20) 



and the magnetic field evolution is governed by the induction equation, derived from the 
Maxw ell equations, and under the considerations listed above (see |Parkeiil2007l : IPandey fc Wardle 
2nn8l . for details). 



dB 

'dt 



= V X 



u X B — ?]J 



Vhall 



jxB + ^(JxB)x 



(21) 



The right hand side of the induction equation has the convective, Ohmic, Hall, and ambipolar 
terms, from left to right respectively. Note that for simplicity we are referring to the Ohmic 
and ambipolar terms as diffusio n terms, but strictly speaking none of them can be cast in 
the form of a diffusion equation ( Parker 1963 1. already used ambipolar diffusion terminology). 
The two new terms (Hall and ambipolar) are implemented in the Bifrost code in the induction 
equation and in the electric field. 

Note that from Equation [211 the Hall and ambipolar terms can be considered as advec- 
tion terms: 



9B 

— = V X [u X B - 77J - uh X B + ua X B] (22) 
at 

where the Hall velocity is Uh = {Vhaii3)/\B\ and the ambipolar velocity is Ua = {Vamb^ x B)/5^ 
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The generalized Ohm's law is implemented in the code using; the same scheme as used 
for the MHD equations, i.e., a 6th order explicit method (IGudiksen et al.l I2OIII ). From 
the expression [22] it is clear that the Hall term and ambipolar diffusivity g ive rise to two 



new constraints on the CFL condition which restrict the timestep interval ( jCourant et al. 



19281 ) {Atn = Ax/uh and At a = Ax/ua). Both velocities are a function of the current 



(V X B), i.e., both CFL conditions are quadratic functions in Ax, and the timestep will 
decrease quadratically with increasing spatial resolution. We note that for the simulation 
with mean magnetic field strength in the photosphere of the order of 100 G, the ambipolar 
and Hall velocities are maximal in the cold regions in the chromosphere with, respectively, 
ua ^ 100 km s~^ and Uh ^ 1 km s~^. As result of this, the CFL criteria are approximately 
AtA ^ 0.3 s and Atn ^ 20 s with Ax ^ 32 km, compared with the strictest CFL condition 
in the simulation of At 3 10~^ s. Therefore, as long as we do not increase the mag- 
netic field and/or the spatial resolution too much, we do not need to change to an implicit 
implementation of our equations. 



2.1.3. The energy equation 

As mentioned above, the energy dissipated by Joule heating is given by Qjouie = E- J. In 
the previous section, the Hall term and ambipolar diffusivity were shown to lead to changes 
in the electric field. These changes need to be taken into account when computing the energy 
due to the dissipation of the magnetic field. Note however, that because the Hall term in the 
electric field is a function of J x B, then (J x B) • J is zero, i.e., the Hall term does not produce 
any energy dissipation at all. The only terms which directly dissipate electromagnetic energy 
by dissipation are by Ohmic and by ambipolar diffusion. In the Bifrost code the former is 
negligible compared to the artificial diffusion needed to stabilize the code at numerically 
resolvable scales and is therefore set to zero. 

In contrast to the artificial resistivity present in the code, the Hall term and ambipolar 
diffusion are calculated as a function of the electron density and, for the latter, of the 
collision frequency between the different species in the solar atmosphere. In order to avoid 
instabilities from rapid heating processes due to the new terms, it is sometimes necessary to 
further limit the time steps (beyond the CFL condition) because the timescales of the energy 
dissipation of the ambipolar diffusion are short. As a result, the energy distribution in the 
chromosphere changes rapidly and the source and sink terms in the energy equation, such 
as radiative processes, need to be updated more often than is the case without ambipolar 
diffusion. 
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2.1.4' Collision frequencies 



The collision frequency between electrons and ions can be found in e.g. iPriestI (119821 ) 
and is given by 



= 3.759 10" VT-^/^ In a 



(23) 



and 



^ = 5.2 10-^1 ^ 



He In A 



(24) 



where In A is the Coulomb logarithm (all in SI units). 



As in iDe Pontieu et al.l ( 120011 ) , we follow three different approxima tions in computin g 



the collision frequency be t ween ions and neutral particles: as descr ibed by 



Osterbrockl (|l961r>: 



De Pontieu fc Haerendell (119981 ) (hereaft er case A) , as de s cribed by I von Steiger fc GeissI ( 119891 ) 
(hereafter case B) and as described by Fontenla et al. ( 19931 ) (hereafter case C), (see Ap- 
pendix Ej). Table [1] lists the 2D simulations for which we investi gate the effects of these 
different methods to calculate Uin- We note that the appendix of ( jPe Pontieu et al.ll200ll ) 
contains two typos: their formula A6 should be divided by 2 to provide the correct expres- 
sion for the collision frequency between neutral hydrogen and protons, and formula A12 
should be replaced by our formula . Our formula provides the correct equation for the col- 
lision frequency between neutral hy drogen and p roton s, according to the recipe derived by 
De Pontieu fc Haerendell ( 119981 ) and lOsterbrockl (ll96ll ) . 



Throughout the paper we will focus on the results of case B since it is more recent, and 
the most extensive. 

In order to calculate the various collision frequencies, the ion and neutral fractions are 
calculated from the Saha-Boltzman equation. The electron density is also computed on the 
basis of LTE, in practice this is done via a table lookup in the Bifrost code. In the pre- 
computed table, the 16 most important atomic species in the solar atmosphere are taken 
into account. Table [2] lists the atomic species, abundances and ionization fraction (Xi). 



2.2. Tests 



One of the main objectives of this work is to study the importance and vahdity of the 
generahzed Ohm's law in a "reahstic" 2.5D simulation of the solar atmosphere. We describe 
three different tests done for the implementation of the generalized Ohm's law, which also 
illustrate the role and importance of each form of diffusivity. For two of the tests, we imposed 
a velocity equal to zero at all times in the full domain. We also run separate tests using only 
the Hall term or ambipolar diffusion. 



2.2.1. ID Hall test 



First, we test that our code correctly includes the Hall term. In this test case, we set 
the velocities and ambipolar diffusion to zero and consider the induction equation in ID 



dBy^ 
dt 

dt 



R ^ 

Vhall^x o 



(25) 
(26) 



For this test, we set constant {Bx = 0, 1121, 2242 G are shown with orange diamond 
symbols, and blue and green lines in Figure [1]). With higher Bx^ the rate at which By and 
Bz change with time increases. However, the total magnetic flux should remain the same 
at all times, since the Hall term cannot convert the magnetic flux into thermal or kinetic 
energy. Note also that the Hall term will give rise to a non-zero B^ (and therefore a non-zero 
Uz in a dynamic simulation) even if the field originally has no component in the z-direction. 
Figure [1] shows By in the top panel and B^ at the bottom panel for four different runs. All 
cases have the same jump in By (black triangle symbols in the top panel) and a constant Hall 
term. In the test shown in red line in Figured] does not have the Hall term and Bx = 2242 G. 
All of these tests are shown at the same instant {t = 20 s). 

The rate of change of By and B^ is as expected, i.e., the case with Bx = 2242 G leads 
to an increase of unsigned total flux of B^ (integrated along the x-axis) that is twice as large 
as the case where 5^ = 1121 G. Moreover, the case Bx = G behaves similarly to the case 
with no Hall term. The magnetic flux is in all cases conserved. This gives us confidence that 
our implementation of the Hall term in the code is satisfactory. 
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2.2.2. ID Ambipolar test 
In ID, the induction equation for By is: 



(27) 



Ot OX \ ^ OX 

Apart from the trivial solution, By = consta nt, it is clear that Equation [271 also permits 
a steady solution of the form of By oc x^/^ (see iBrandenburg fc Zweibell 11994 for details). 



In this expression we should keep in mind that the code includes numerical diffusivity in 
addition to ambipolar diffusion. We consider the evolution of an initially sinusoidal profile 
of By. This profile evolves, and strong gradients become stronger approaching the form 
By OC x^/^ as time progresses. Figure [2] shows the initial condition of By (solid line) and at 
t = 50 s (dashed line) which is close to the steady solution. Observe that where the gradient 
of By is large By closely follows the expression By oc x^/^ (dash-dotted line). 

Ambipolar diffusion converts magnetic energy into thermal energy as discussed above. 
In this ID test, we turn off the heating from artificial diffusion and only allow heating from 
the ambipolar diffusivity. Such heating in this simple simulation must follow the expression 



^ = VambJ'.Bl (28) 

Figure [3] shows the energy profile with x of this test (black diamonds) at t = 2.1 s. 
Calculating the right hand side of this expression using the sinusoidal shape of By from the 
initial condition, then deriving J^, we calculate the energy to be 



e = einit + TJambJzByAt, (29) 

where At is the time increment. This relationship is shown with the red line in Figure [3l 
The black diamonds overlaps the red line as would be expected. This indicates we have 
correctly implemented ambipolar diffusion in the code. 



2.2.3. Collision frequencies test: VAL-C model 

In order to test whether the absolute values of the diffusion terms are calculated cor- 
rectly, we use three different sources for the neutral-ion collision frequency (z/^^, see Sec- 
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tiQn l2.1.41) . This also allows us to study the uncertainties involved i n the various forniulas fo r 
the collision frequency, as already studied for ID-static models by iDe Pontieu et al.l ( 120011 ). 
We test our implementation in t he Bifrost code by us ing the densities and temperatures from 
the VAL-C atmospheric model (IVernazza et al.lll98ll ) which allows us to compare our results 
with those found in the published literature. Indeed, we correctly obtain the neutral-ion 
collision fre quencies as a func t ion o f height as can be seen by comparing our Figure H] with 
Figure 2 in iDe Pontieu et al.l ( 120011 ). The large dip in the collision frequency at 0.5 Mm is 
due to low number of ions (mostly non- hydrogen species) in this region.. 



2.3. Initial and boundary conditions 

Let us now consider the importance and validity of the generalized Ohm's law in a "re- 
alistic" 2.5D simulation of the solar atmosphere. The 2.5D computational domain stretches 
from the upper convection zone to the lower corona and is evaluated on a non-uniform grid 
of 512 X 325 points spanning 16 x 16 Mm^. The frame of reference for the model is chosen 
so that X is the horizontal direction and z is the vertical direction (Figure [5j). The grid is 
non-uniform in the vertical z-direction to ensure that the vertical resolution is good enough 
to resolve the photosphere and the transition region with a grid spacing of 28 km, while 
becoming larger at coronal heights where gradients are smaller. 

We run two different initial conditions with different values for the unsigned magnetic 
field strength but with similar field configurations (Figure [5j). By is originally set to zero. 
The initial model starts with a magnetic field that is inclined some 5 degrees with respect 
to the vertical axis and the two different setups for the unsigned field strengths in the 
photosphere are 0.25 G, the other 90 G. These two initial conditions are run for the three 
different formulas that were mentioned above to calculate the collision frequency z/^^. The 
simulations with the initially weak magnetic field using cases A, B or C for the neutral-ion 
collision frequency are labeled WA, WB or WC, while the strong field simulations using cases 
A, B, or C for the neutral-ion collision frequency are labeled SA, SB, or SC, respectively. 
Tabled] list the different simulations. 

In the following we will refer in our analysis to simulations WB and SB, unless otherwise 
noted. 
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Results 



The simulations presented include the dynamic processes (including radiative losses) of 
the photosphere and chromosphere and a self-maintained chromosphere and corona. This is a 
very different type of model as compared to semi-empirical models such as the VAL-C model, 
or previous simulations investigating the effects of partial ionization which had a simplified 
treatment of the energy balance (and ionization degree). It is thus of significant interest to 
determine how dynamic atmospheres such as from our simulations impact the importance 
of the ambipolar diffusion and the Hall terms (also assuming different approximations to 
the collision frequency), and to compare the results with those from the models based on a 
VAL-C type atmosphere. 

The basic structure of our modeled chromospheres are shown in Figure [51 A full de- 
scription of their properties fall outside the scope of this paper but we will mention the most 
important as concerns ambipolar diffusion and the Hall term: the basic thermodynamic state 
of the chromosphere is maintained by the continual injection of acoustic shocks from the pho- 
tosphere. These perturbations are due to the chaotic generation of waves in the convection 
zone, of which waves with periods of order 3 minutes will propagate and steepen in th e 
chromosphere, as is well known and as extensively studied by lCarlsson fc Stein (ll992L Il994r). 



The propagation of waves will be modified in the pre s ence of a magnetic fie l d (jBogdan et al 



2011 



20031 : iDe Pontieu et al.l l2004l : iHeggland et all 120071 : Mcintosh et all l201ll : iHeggland et al. 



, among others) 



but will nevertheless steepen and form strong shocks, with high temperature s in the 
shock fronts and very low temperatures in the regions behind (ILeenaarts et al.ll201ll ). These 
"cold chromospheric bubbles" can be seen in both panels of Figure [5] which show temper- 
atures as low as 2 000 K or lower. In the strong field case the Lorentz force is clearly 
important, pushing the corona upwards and allowing cool material to exist at great heights, 
much higher, up to 5 Mm above the photosphere, than that found in semi-empirical models 
where the maximum chromospheric height is found to be of order 2-2.5 Mm. The distribu- 
tion of density and temperatu re with height in dynam ical "realistic" simulations is discussed 
in much greater detail in e.g. ILeenaarts et al.l ( 120111 ). 



3.1. Collision frequencies and diffusivities 

As mentioned, most studies of the effects of ion-neutral collisions in the chromosphere 
have been based, in some form, on semi-empirical models (VAL-C or FAL-C models as shown 
in Figure H]). However, the chromosphere and transition region are clearly highly dynamic. 
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and it is of great importance to know the effects of the neutral-ion interactions in such 
dynamic atmospheres. First of aU, we are interested in studying the relative importance of 
the different diffusivities in the chromosphere and transition region. Figures [6] and [71 show the 
Ohmic diffusion, artificial diffusion, Hall term, and ambipolar diffusion from top to bottom 
and left to right for the simulations WB and SB, respectively. 

On comparing the different diffusivities, we find that in the entire chromosphere, the 
Hall term is on average two orders of magnitude larger than Ohmic diffusion. This is true 
for both simulations WB and SB and is perhaps more easily seen by considering the ratios of 
the diffusivities plotted in Figure [81 Ambipolar diffusion is roughly four orders of magnitude 
larger in the weak field (WB) case and fully six orders of magnitude larger than Ohmic 
diffusion for the strong field SB case. Ambipolar diffusion is considerably larger for SB than 
for WB because ambipolar diffusion depends quadratically on the magnetic field strength. 
Note that while Ohmic diffusion has a significant magnitude throughout the atmosphere, 
ambipolar diffusion is important only in the chromosphere. 

Numerical simulations must include some form of artificial diffusion in order to com- 
pensate for the fact that they do not have infinite spatial resolution. In the Bifrost code 
this is done through a so called hyper-diffusivity which in practice means that the diffusion 
coefficient is increased in regions that require high diffusivity, i.e., where gradients are large. 
The magnitude of this artificial diffusion is set by the spatial resolution. To some degree this 
behavior is similar to Ohmic diffusion, but there are also significant differences. Simulations 
run at the highest possible spatial resolution cannot even come close to the diffusion values 
found in Ohmic diffusion. This is a well known problem for numerical simulations of the 
solar atmosphere. 

For the simulations reported here, we find an artificial diffusivity in the chromosphere 
that is three orders of magnitude larger than the Ohmic diffusivity and that is up to five 
orders of magnitude larger than the Ohmic diffusivity in the corona. By design, the artificial 
diffusion is largest where the shocks and other high gradient phenomena are located. More- 
over, since the grid is non-uniform and the grid resolution decreases with height, artificial 
diffusion will on average be larger in the corona than in the lower layers. In contrast, the 
Ohmic diffusion is largest in the upper photosphere and chromosphere. As result of these 
differences between artificial diffusion and Ohmic diffusion, the magnetic Reynolds number 
on the Sun and in the simulations is completely different: the magnetic Reynolds number is 
several orders of magnitude higher in the solar atmosphere than even the highest resolution 
simulations. Since Ohmic diffusion is negligible compared to artificial diffusion we do not 
include its effects either in the induction equation nor in the energy equation. In a similar 
manner, as result of the low resolution of these simulations, the artificial diffusion will mask 
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the Hall term, but here we are interested in describing the stratification of the Hall term and 
ambipolar diffusion in self-consistent magneto-convection simulations. 

One of the most interesting results of our calculations is that, on the other hand, am- 
bipolar diffusion is of the same order or, in some regions, even larger than the artificial 
diffusion in the chromosphere. This is a perhaps surprising, but crucial property of the chro- 
mosphere. It allows us to use our numerical simulations to study the effects of ambipolar 
diffusion while using the correct physical magnitudes of the coefficients. As a result, the 
chromosphere may be the only region where simulations are close to reality once all the 
physics are included in the code, despite the necessarily limited resources of todays com- 
puting technology. This has an impact beyond the chromosphere, since it directly affects 
discussions on whether these self-consistent magneto-convective simulations provide a real- 



i stic d river and boundary to the corona. For example, recent simulations bv lHansteen et al. 



( I2OIOI ) suggest a preponderance of heating in the lower atmosphere (first few Mm above 



the photosphere), implying th at much of the coronal heating occurs towards the footpoints 



( iMartmez-Sykora et al.ll201ll ). The large ambipolar dissipation we find here suggests that 
such simulations (which only include artificial resistivity) are actually much more realistic 
than previously thought, including the predictions of heating low down. 

The Ohmic diffusivity. Hall term, and ambipolar diffusivity depend on the electron den- 
sity, while the Ohmic and ambipolar diffusivites also depend on collision frequencies which 
are shown in Figures [9] and [10] for the weak field WB and strong field SB simulations (see 
Appendix Ej). (Note that the Ohmic diffusivity is proportional to the collision frequency, 
while ambipolar diffusivity is inversely proportional to the collision frequency.) On aver- 
age, the collision frequency between electrons and ions is larger than both the ion-neutral, 
electron-neutral, and neutral-ion collision frequencies in the chromosphere, in the WB simu- 
lation one order of magnitude larger and in the SB simulation two orders of magnitude. This 
difference in collision frequencies between the WB and SB simulations is mainly because the 
chromosphere is hotter in the SB simulation as a result of ambipolar heating. 

The electron-ion collision frequency also shows strong variation throughout the chromo- 
sphere, by almost 5 orders of magnitude in both simulations. This variation is due to the 
electron density variation in the chromosphere (see second row in Figured!] and Equation [231). 
As a result, the electron-ion collision frequency is lower inside the cold chromospheric bubbles 
than in the shock fronts. In the chromosphere, the electron-neutral and ion-neutral collision 
frequencies are similar in magnitude. However, the ion-neutral collision frequency shows a 
stronger variation in space in the middle chromosphere than the electron-neutral collision 
frequency. This is especially true in the cold chromospheric bubbles, where the ion-neutral 
collision frequency is significantly lower. What causes these differences? First, we note that 
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Pn shows less variation in horizontal cuts in the lower chromosphere than pi because the re- 
gion is mostly dominated by neutrals. As result of this, the neutral density is almost similar 
to the total density. In the cold bubbles, hydrogen is mostly neutral, and the only ions are 
provided by the heavier metals. While both the electron-neutral and ion-neutral collision 
frequency are dependent on the neutral density (which does not vary much in the lower 
chromosphere), the dominance of metals in providing ions implies that the average mass per 
ion increases significantly in the cold bubbles (compared to the rest of the chromosphere). 
The associated drop of average thermal speed (for the heavy ions compared to protons) is 
the reason for the sharp drop in ion-neutral collision frequency in the bubbles (compared 
to the rest of the chromosphere). The neutral-ion collision frequency is even lower than 
the ion-neutral collision frequency in the bubbles, because there are so few ions available to 
collide with (bottom panels in Figure [H]) . 

We now consider which parameters are responsible for the changes in diflFusivities through- 
out the solar atmosphere. In both simulations (WB and SB), the strongest Ohmic diffusivity 
is concentrated in the lower-middle chromosphere while it is weaker in the corona and con- 
vection zone. In the chromosphere, the Ohmic diflFusivity varies over a range of almost four 
orders of magnitude. This variation in the chromosphere is due to the strong variation of 
the electron density and collision frequency of electrons with neutrals and ions (Figures [9]- 
[TT]) . Ohmic diffusion is large in the expanding cool bubbles and low where temperatures 
are higher. This is because the Ohmic diflFusion variations are dominated by the variations 
in electron density, which is very low in the cool bubbles, and large in shock fronts. The 
collision frequency of electrons with ions and neutrals does not drop as precipitously in the 
cold bubbles since there are plenty of neutrals to collide with in these bubbles. 

The Hall term is largest in the lower-middle chromosphere and in the corona (Figure [8]). 
This is because it is inversely proportional to the electron density which is small in both 
regions. We see that for both simulations, the Hall term is larger than the Ohmic diffusivity 
in the chromosphere and corona, but not in the photosphere nor in the convection zone. In 
the cooler regions of the chromosphere, the Hall term is relatively even higher than in the 
shock fronts, and up to three orders of magnitude greater than the Ohmic diffusivity. Such 
differences are a bit larger in the WB simulation, since electron density is smaller in the 
cold chromospheric bubbles in the weak field model. This difference in the electron density 
between WB and SB is because the cold bubbles have cooler temperatures in WB simulation 
than in the simulation SB. In the intergranular lanes in the photosphere, the Hall term is the 
most important diffusion term after the artificial diffusion. Therefore, since the Hall term is 
proportional to the strength of the magnetic field, this term may be important to consider 
in magnet oconvective simulations that include strong magnetic fields. 
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Ambipolar diffusion is important in the region from the upper-photosphere to the upper- 
chromosphere. In the photosphere, ambipolar diffusion shows some importance in intergran- 
ular lanes which have strong concentrations of magnetic field (Figure [H]). Therefore, the 
strong field in the SB simulation shows considerably more diffusivity in intergranules with 
high magnetic flux concentrations than in the weak field WB simulation. In the chromo- 
sphere, ambipolar diffusion dominates almost everywhere except for in the lower chromo- 
sphere in shock fronts. The largest difference between ambipolar and Ohmic diffusion is 
located in the cold chromospheric bubbles and near the upper-chromosphere/lower tran- 
sition region. Note that for the SB simulation the ratio between ambipolar and Ohmic 
diffusion is almost four orders of magnitude larger than in the WB simulation due to the 
quadratic dependence of the ambipolar diffusivity on the magnetic field strength. The am- 
bipolar diffusivity is large in the cold bubbles since the ion density and the ion-neutral 
collision frequencies are low, but mainly because the ion density is extremely low (5 orders 
of magnitude lower than in the chromospheric shock fronts). In the upper chromosphere, 
the ambipolar diffusivity becomes relatively strong due to low densities — which lead to 
low ion-neutral collision frequencies — but only in those regions where the magnetic field 
strength is high. In the cold bubbles the ion density is low because of the adiabatic ex- 
pansion and cooling, whereas in the upper chromosphere, it is because the density drops 
by 2 3 orders of magnitude compared to the lower chromosphere. As shown in Figure [HI 
the ambipolar diffusivity depends strongly on the ion and neutral density and thus on the 
ionization state of the chromospheric plasma. However, it is well known that in the middle 
and upper chromosphere the ionization and r ecombination rates are fairly slow for hydrogen 
which will not be in ionizational equilibrium ( jCarlsson &: Steinll2002l ). This suggests that, in 
order to treat ambipolar diffusion realist ically, it is necessary t o solve the full time dependent 
rate equations for hydrogen ionization ( ILeenaarts et al.l 120071 ). 



3.1.1. Comparison with VAL-C model 

The VAL-C model does not provide a good description of the strong temporal and 
spatial variations found in the physical variables of the chromosphere. In the section above, 
we have seen that the different collision frequencies and diffusivities show spatial variations of 
several orders of magnitude at the same height in the chromosphere. The lower- chromosphere 
changes rapidly due to the shock fronts; these lead to changes in the thermodynamic structure 
of the lower chromosphere on times-scales shorter than a minute. As a result of this, cold 
chromospheric bubbles appear and disappear in minutes. Due to the ambipolar diffusion, 
plasma is heated in the cold bubbles on timescales shorter than those characterizing the shock 
front. As a result of the spatial and temporal variations, the neutral-ion collision frequency 
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varies by almost eight orders of magnitude in the chromosphere in the 2D simulation, whereas 
the VAL-C model has a unique value for the collision frequency at every height (Figure [T2l). 
In the cold chromospheric bubbles, the collision frequency drops to considerably lower values 
than those found in the VAL-C model. This is a result of the low ion number density in 
these areas, which are overestimated in the VAL-C model. 

We use the maximum, minimum and median magnetic field of the 2D models (SB and 
WB) as a function of height in order to calculate the range of ambipolar diffusivities in 
the semi-empirical VAL-C model. The ambipolar diffusion has a very wide range of values, 
8 orders of magnitude in simulation SB and 11 orders of magnitude for simulation WB. 
These variations are almost 6 or 8 orders of magnitude larger than in the VAL-C model. 
The ambipolar diffusivity in the 2D simulations is much higher in most of the chromosphere 
compared to what is found in the VAL-C model. The reason for the large difference of the 
neutral-ion collision frequency and the ambipolar diffusivity in the VAL-C model (compared 
to the 2D model) is because the VAL-C model does not capture the thermodynamics of the 
cold chromospheric bubbles where the neutral-ion collision frequency drops precipitously. 

These large differences in both the neutral-ion collision frequency and ambipolar diffu- 
sivities, found between the VAL-C model and our simulation should lead to a re-examination 
of previous results related to the generalized Ohm's law (see references) using semi-empirical 
models to define the density and temperature structure. We also reit erate the importance o f 



taking into account the likely dynamic state of hydrogen ionization ( ILeenaarts et al.l 120071 ). 



3.1.2. Other methods to calculate collision frequencies 

We considered three different methods to calculate (Section 12.1.41) . and thus, the 
ambipolar diffusivity. Do the different methods give similar values of the collision frequency 
and/or diffusion for the different models? We note that the evolution of the simulations 
using the different methods to calculate the collision frequency (WA, WB and WC, and SA, 
SB, and SC) diverge within a few minutes. We therefore integrate the properties of the 
models in time in order to study the different values of the collision frequency and ambipolar 
diffusion, and proceed as follows. In Figure [13] we show the joint probability distribution 
function (JPDF) of the temperature versus the ion-neutral collision frequency (top) and the 
ambipolar diffusivity (bottom) for the simulations labeled WB (left) and WC (right). In 
Figure [14] we show the cases SB and SC in a similar manner. We have integrated over 4 
minutes. Note that the variation of the ^-axes are logarithmic and cover more than 10 orders 
of magnitude. We mostly focus on cases B and C since those are the most recent, and include 
more advanced calculations of the collision frequencies. 
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The values for and rjamb differ in range and mean values for the different cases in each 
simulation. These differences are significant in certain temperature ranges. The differences 
between the different methods to calculate the ion-neutral collision frequency are similar for 
both atmospheres (weak and strong magnetic field strength). For instance, at log(T) ^ 3.7 
(5000 K), case B shows ion-neutral collision frequencies that are a factor two larger than for 
case C. At temperatures larger than log(T) ^ 3.7, the collision frequency for case B is almost 
two orders of magnitude smaller than in case C. As a result, in certain temperature ranges, 
the largest values of the ambipolar diffusion for case B are almost 2 orders of magnitude 
larger than for case C. For temperatures lower than log(T) ^ 3.6 (4000 K), the median 
collision frequency as a function of temperature is roughly similar between cases B and C, 
but not the distribution, as can be seen: case B reaches collision frequencies smaller than 
case C. 

At temperatures larger than log(T) ^ 3.8 (6300 K), the collision frequencies for case C 
are one order of magnitude larger than for case B. As a result, the median of the ambipolar 
diffusion for cases WC and SC is one order of magnitude smaller than for WB and SB. 

In order to have a better impression where in the atmosphere the ion-neutral collision 
frequency and ambipolar diffusion differ between the different cases, we take the same atmo- 
spheric model (simulation WB or SB at t = 2500 s) and calculate from these two models the 
collision frequency and ambipolar diffusivity using the different methods (Figures [T5lfT6l) . It 
is interesting and important to see that at the precise location where the ambipolar diffusion 
is really high (in cold chromospheric bubbles and in the upper chromosphere), the different 
methods differ most. In the cold chromospheric bubbles for both atmospheres (weak and 
strong magnetic field), the collision frequency using case B is almost 4 times smaller than 
case A, but similar to case C. As result, the ambipolar diffusivity is more than 2 times smaller 
using the method of case B than for case A. In the upper chromosphere or shock fronts, the 
collision frequency using case B is almost two times larger than case A and slightly larger 
than with case C. As a result, the ambipolar diffusivity using case B is more than 3 times 
smaller than case A, and almost 3 times smaller than case C. However, in the proximity of 
the transition region, the collision frequency for case B is a bit smaller than case A, and 
more than 10 times smaller than case C. As result of this, the ambipolar diffusivity using 
case B is a bit larger than case A, and more than 10 times larger than case C. 

These large differences between each method are due to the different temperature de- 
pendences (see Appendix Ej). As mentioned, these differences lead to rapidly diverging ther- 
modynamic evolution in the various models. Thus, it is important to take into account 
this uncertainty in calculating the collision frequencies when the generalized Ohm's law is 
modeled. 
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3.2. Approximations to the generalized Ohm's law 

The generalized Ohm's law is based on several approximations and considerations. In 
this section we describe where these approximations fail and the implications of this failure. 
We employ the atmospheres of the WB and SB simulations in this discussion. 



3.2.1. Approximations in the momentum equation 

Let us establish and validate the different assumptions underlying the generalized Ohm's 
law as implemented in the code, and see if they are fulfilled in the fully dynamic self- 
consistent simulations. One of the first consideration is that the ion density dominates over 
the electron density {pi/p » Pe/p)- Everywhere in the atmosphere, the values of ion and 
electron densities remain within the range that fulfill pi/p » p^/p so that electron inertia 
can be neglected. 

In order to neglect the effects of drift momentum in the momentum equation, the drift 
momentum has to be s maller than the fast momentum {p-s/vl + c^, see Equation [6] and 



Pandey fc Wardld (120081 )). This approximation is fulfilled in most of the atmosphere under 



both strong and weak field conditions. The only exception is in the weak field atmosphere, 
where some low density areas just below the transition region show a ratio of order 0.1-1 
(see Figure [TTj). This is because the ion-neutral collision frequency drops significantly there, 
so that the drift between ions and neutrals becomes rather large. As a result, in these small 
regions the plasma becomes decoupled from the neutrals, and it may be necessary to add 
the drift momentum to the momentum equation, and/or solve the MHD equations using 
multiple fiuids. In the weak field atmosphere, a few of the cold, expanding bubbles show 
ratios of order 0.1, so that the fast momentum does stay in excess of the drift momentum. 
This suggests that the region below the transition region is the only one of concern for this 
particular condition. 



3.2.2. Approximations in the induction equation 



To allow the removal of the time dependence of the drift momentum equation ( jPandey fc Wardle 



20081 ) . the electron density times the collision frequency of electrons with neutrals has 



to be smaller than the ion density times the collision frequency of ions with neutrals ( 
Pe^en << Pi^in) • This approximation is fulfilled in most of the atmosphere with the ex- 
ception of some areas in the upper photosphere and in the cold chromospheric bubbles 
(Figure [T8l) . In the cold bubbles, the electron-neutral collision frequency is almost similar 
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to the ion-neutral collision frequency. In the upper photosphere, and the collision frequency 
of electrons with neutrals is relatively large so Pe^en « Pi^in is not fulfilled. Therefore, in 
these regions, the proper way to solve the ambipolar term in the induction equation is by cal 



culating the drift velocity using the fully time dependent equation of Ud ( jPandey fc Wardle 



20081). 



Some of the approximations used in deriving the equations require that the dynamical 
frequency remains smaller than the frequencies shown in Figure [191 The typical timescales 
on which the simulated atmosphere evolves is of order 10s or longer, i.e., a dynamic frequency 
of ^ 0.5 Hz or lower: if the frequencies shown in Figure [19] are higher than ^ 0.5 Hz, the 
assumptions underlying the generalized Ohm's law are fulfilled. 

The first assumption is that the time derivative of the drift velocity can be neglected. 
Following Equation [181 this can be done only if the dynamical frequency is smaller than 
the frequency {p/ Pi)i^ni shown in the top panels of Figure [T9l The latter frequency is very 
high in the upper photosphere and the chromosphere, and stay well above the dynamical 
frequency of our simulations {'^ 0.5 Hz). Only in the vicinity of the transition region does 
{p/ Pi)yni become small enough that it is of the same order as the dynamical frequency of 
the simulations. As a result, we may need to take into account the derivative terms shown 
in Equation [T7] only in this small region in the vicinity of the transition region. 

A second assumption i s that the dyadic produ ct of the drift velocity in the momentum 



equation can be neglected ( iPandey fc Wardldl2008l ). This term can only be neglected if the 



dynamic frequency stays well below the frequency defined in Equation [7] and shown in the 
middle panels in Figure [T9l We find that in the cold chromospheric bubbles (in the weak field 
case) and in the upper chromosphere (in both weak and strong field cases), this assumption 
sometimes fails. In these regions, we may thus need to take into account the momentum 
drift term in the momentum equation. 

A final assumption is that the terms of the form VP x B in the induction equation can 
be neglected. This can only be done when the dynamic frequency stays below the frequency 
defined in Equation [19] and shown in the bottom panels in Figure [191 This bound for the 
dynamical frequency strongly depends on magnetic field strength of the model. We find that 
for the weakly magnetic atmosphere case (WB) this limit is low, and the assumption fails 
in the upper chromosphere and cold bubbles. In the strongly magnetic atmosphere (SB) the 
assumptions only fails in the upper part of the chromosphere. 

In summary, for the weak field atmosphere we cannot neglect the time derivative and 
dyadic product of the drift velocity, and the VP x B terms (in the momentum and induction 
equation respectively) in the cold bubbles and just below the transition region. In all other 
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regions in the weak field atmosphere, the assumptions underlying the generalized Ohm's law 
are fulfilled. 

For the strong field atmosphere, the generalized Ohm's law works well in most of the 
chromosphere, except in the region just below the transition region where the time derivative 
and dyadic product of the drift velocity and the VP x B terms cannot be neglected. 
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4. Discussion and Conclusions 

We have implemented the partial ionization effects in the Bifrost code in the form of the 
Hall term and ambipolar diffusion. The code has been tested and verified with different tests 
that are presented in this paper. The code allows the simulation of the solar atmosphere, 
from the upper convection zone to the lower corona, with a magnetoconvective photosphere, 
and a fully-dynamic and self-maintained chromosphere and corona. We studied the different 
diffusivities in two different models, one is weakly magnetic, and the other is rather strongly 
magnetic. The magnetic field strength of the latter model is similar to that found in the 
quiet sun, including the network. 

In short, the Ohmic diffusion is roughly three orders of magnitude smaller than the 
Hall term in the chromosphere, and the latter is three orders of magnitude smaller than the 
artificial diffusion. Unlike Ohmic diffusion, the Hall term depends on the magnetic field, 
as does ambipolar diffusion which is strongly dependent on the magnetic field strength. As 
a result of this, the ambipolar diffusivity is clearly different for the two models; in regions 
with large ambipolar diffusivity we find it is of the same order as the artificial diffusion 
in the chromosphere for the weakly magnetic model (WB), and more than one order of 
magnitude larger than the artificial diffusivity for the strongly magnetic model (SB). The 
fact that the artificial diffusivity is actually smaller than the ambipolar diffusivity under 
many chromospheric conditions has some very important consequences. It means that these 
simulations are capable of providing a surprisingly realistic view of the consequences of 
the ambipolar diffusion in the chromosphere and corona. This has an impact beyond the 
chromosphere, since it directly affects discussions on whether these self-consistent magneto- 
convective simulations provide a realistic driver and boundary to the corona. These results 
will be described in detail in a follow up paper. 

Another important result is that both, the Hall term and ambipolar diffusivity, vary 
by several orders of magnitude in the chromosphere as result of the time varying dynamics 
and the strong variations in temperature, electron, ion and neutral density, and magnetic 
field strength in this region. This strong variation is not taken into account in any of 
the previous studies which use either ID semi-empirical VAL-C type models, or lack more 
sophisticated approaches to the radiation, ionization and energy balance. The largest values 
of the ambipolar diffusivity are located in the cold chromospheric bubbles that have low 
temperatures due to strong adiabatic expansion, and in the upper chromosphere because 
the neutral-ion collision frequency is small. Howe ver, the ambipolar diff usion is strongly 



dependent on the ionization degree, and as shown by lLeenaarts et al.l ( 120071 ). time dependent 
hydrogen will change the ratio between neutrals and ions compared to LTE conditions. 
The Bifrost code can treat the time-dependent ionization of hydrogen and we plan to run 
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new simulations taking into account both the generahzed Ohm's law and time-dependent 
hydrogen ionization. 

We have compared different methods to calculate the collision frequency between neu- 
trals and ions. Both the ion-neutral collision frequency and ambipolar diffusivity differ 
considerably as a function of the method used to calculate this collision frequency. Since 
ambipolar diffusion has a significant impact on the thermodynamic evolution of these models, 
the simulations rapidly diverge. When comparing each method we find the largest differences 
are located in regions where the ambipolar diffusivity is large: in the cold chromospheric bub- 
bles and in the upper chromosphere in the vicinity of the transition region. These differences 
bring a new uncertainty to the results (Section 13.1.21) . and highlight the need for a detailed 
consideration of the relevant collisional processes in the chromosphere. 

Finally, we investigated t he different approx i matio ns underlying the generalized Ohm's 



law as described in detail bv IPandey fc Wardld ( 120081 ) . In both models, most of the sim- 



plifications are applicable with some exceptions. In the upper-chromosphere the collision 
frequency is too low, as a consequence, the velocity drift can be large. Therefore, we may 
need to define the velocity drift and add an extra term in the momentum equation related 
to the momentum drift between ions and neutrals. In the upper photosphere, and in cold 
chromospheric bubbles the ambipolar term in the induction equation may need to be cal- 
culated using the drift velocity. M oreover, the drift velocity should be calculated using the 



time dependent form (as shown in IPandey fc Wardld l2008l ) . This is necessary because the 



ion density and the ion-neutrals collision frequency drop in these cold areas as opposed to 
the electron density and the electron-neutral collision frequency. 
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A. Collision frequencies 



In order to calculate the collision frequency between ions and neutral particles we use 



three diffe rent approx i matio ns (following the approach by iDe Pontieu et al.l 120011): one de - 
scribed by lOsterbrockl (jl96ll ) ( hereafter case A) , one described bv lvon Steiger fc GeissI (119891 ) 
(hereafter case B) and one by iFontenla et al.l (119931 ) (hereafter case C). 

As a first approach (case A), we take the formulas from Qsterbroc^ ( 1961 ) and De Pontieu fc Haerendel 



(119981 ). where the collision frequency between neutral hydrogen and protons {uhp) are given 
by: 



1 l8kT 



-«P = 510-'%/5^^np (Al) 



rriH is the hydrogen atom mass, and Up is the proton number density. Note that lDe Pontieu et al 



( I2OOII ) had a typo with a factor of 2. The collision frequency {i^Hm) of neutral hydrogen with 



an ionized metal is defined as: 



^Hm = 8 lQ-'\l^^^l^n^ (A2) 
V + 1 V vrm/f 

where and are the atomic mass number of metals ions and the number density of 
metals ions of type m, respectively. The collisions between neutral helium and ions is given 
by: 



This preprint was prepared with the A AS lATgX macros v5.2. 
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UHep = 4 10-^\ - J Up 

V 5 V vrmn 



PHep = 4 10 



-20 



rrim + 1 V TTTTlH 



8A:r 



(A3) 
(A4) 



For the second approach (Case B), foUowing lDe Pontieu et ahl (l200l[ ). lvon Steiger fc Geiss 



( 1l989r i describe the colhsion rate as follows: 



118W^ ( 1-0.125 log^l ^ 



UHm = 21.05 



A~+T '"10^ 



(A5) 
(A6) 



For the helium-proton and hehum-metal coUision frequency we foUow iGeiss fc Buergi 



ep 



Zri 



5.84 



-Zjrr 



(A7) 
(A8) 



where is the ionization weig ht and we conside r ed th at the ions have only one ionization 
state, i.e., = 1. Note that iDe Pontieu et al.l ( 120011 ) have a typo where the expression 
for Uffep is missing the square root symbol for the temperature and the constant 2.2 is also 
different. 



Finally, we find the collision frequencies for the Case C in the appendix of lFontenla et al. 



(119931). 



Using these collision frequencies (Eq. lAlHASp . the collision frequency of neutral Hydro- 
gen with all ions is given by: 



T^Hi = ^Hp + ^HC + ^HN + ^HO + ^HNe + ^HNa + ^HMg + ^HAl (A9) 
+ l^HSi + ^HS + ^HK + T^HCa + T^HCr + T^HFe + T^HNi (AlO) 
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and similarly for the collision frequency of neutral Helium with all ions. Finally, the average 
neutral-ion collision frequency is given by 

PH . pHe /A11\ 

l^ni = ^Hi H ^Hei (All) 

Pn Pn 

Note that in the main text we often use z/^^, which can be derived from using 
momentum conservation {pji^jk — Pk^kj)- 
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Fig. 1. — By (top panel) and (bottom panel) as a function of x are shown for the different 
ID simulations with constant Hall term at time t = 20 s. The initial condition is the same 
for all simulations (shown with black triangles). The runs have different constant Bx values: 
Bx = G (orange diamonds), = 1121 G (blue line), B^ = 2242 G with the HaU term 
(green line), and Bx = 2242 G without the Hall term (red line). Note that the orange 
diamonds, red line, and black triangles overlap. 
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Fig. 2. — From the ambipolar test. By is shown as a function of x at t=50 s. The initial 
condition is shown in solid line. The dashed line shows By dXt = 50 s. The dash dotted line 
shows a function proportional to x^/^ which is what would be expected from the analytical 
considerations. 
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80 




Fig. 3. — Test of ambipolar diffusion on the energy balance. Energy is shown as a function of 
X at t=2.1 s. The energy from the model is shown with the black diamonds and the energy 
extracted from Equation [29] is shown with the red line. Note that the red line is overlapping 
with the black diamonds. 
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Fig. 4. — Neut r al-ion collision frequency as a function of height for the quiet sun model of 
Vernazza et al.l ( ll98ll ). using different formulas for z/^^: the dotted line is case A, solid line 
is case B and dashed line is case C. 
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Fig. 5. — 2D snapshots of the two initial 2.5D MHD models. The initial conditions with 
weak (simulations labeled WA, WB, and WC) and strong magnetic field (SA, SB, and SC) 
are shown respectively in the left and right panel. The color scale shows the temperature in 
logarithmic scale and the magnetic field is shown with white lines. 
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Fig. 6. — Comparison of the different diffusivity terms for the simulation WB at t = 500 s. 
Vohrm Vart) Vhaih Vamb ^re shown from top to bottom and left to right respectively in 
logarithmic scale. Note that more than 9 orders of magnitude are shown. 




Fig. 7. — Comparison of the different diffusivities for the simulation SB at t = 500 s. The 
layout is the same as in Figure [6l 
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Fig. 8. — Ratio between the Hall term and Ohmic diffusion (left panels) and ambipolar and 
Ohmic diffusion (right panels) for the SB simulation (top panels) and WB (bottom panels) 
at t = 500 s. 
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Fig. 9. — Comparison of the different collision frequencies for the WB simulation at t = 500 s. 
Uei^ Pen) ^ud u^i are shown from top to bottom and left to right respectively in logarithmic 
scale. 
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Fig. 10. — Comparison of the different collision frequencies for the SB simulation at t = 500 s. 
The layout is the same as in Figure M 
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Fig. 11. — The ambipolar diffusion, electron density, absolute value of the magnetic field, 
ratio between neutral and total density and ion density are shown from top to bottom for 
the WB simulation (left panels) and SB (right panels) at t = 500 s. 
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Fig. 12. — Minimum (dashed line), median (solid line) and maximum (dashed line) of (top 
panels) and rjamb (bottom panels) as function of height are shown for the simulation labeled 
WB (black in left panels), SB (black in right panels) and for the VAL-C atmosphere (red). 
The VAL-C ambipolar diffusion is calculated taking into account the maximum, minimum 
and median magnetic field of the 2D models as a function of height. The minimum, median 
and maximum are calculated in horizontal planes for the instant t = 500 s. 
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Fig. 13. — Joint probability distribution function (JPDF) of the temperature against of the 
ion-neutral collision frequency (top) and ambipolar diffusivity (bottom) for the simulations 
labeled WB (left) and WC (right). JPDF is calculated integrated over 220 s above the 
photosphere. The colorbar is in logarithmic scale. The median as a function of temperature 
for the WB and WC cases are shown in solid, and dashed lines respectively. 
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Fig. 14. — The layout is the same as Figure [131 However, the simulations are SB (left) and 
SC (right). 
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Fig. 15. — Ratio of the ion-neutral collision frequency (left panels) and ambipolar diffusivity 
(right panels) between case A to case B (top panels) and case C to case B (bottom panels), 
for the atmosphere with weak magnetic field strength. The white contours show where these 
ratios are equal to one. Note that the color scheme is in a logarithmic scale. We used the 
same atmospheric model for all three cases, i.e., before the simulations diverge with time, but 
then calculated the collision frequency and ambipolar diflFusion using the diflFerent formulas 
of each case. 
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Fig. 17. — The drift momentum has to be smaUer than the fast momentum (Equation [6]). 
The ratio between both terms, i.e., p^p^u^ and p^{vl^ + c^) is shown for the simulations 
labeled WB (top panel) and SB (bottom panel) at t = 500 s. The colorbar is in logarithmic 
scale and it is the same for both panels. 

Table 1. Simulation description 



Name Collision frequency Min/Mean/Max \B\ [G] 



WA Case A 0.003/0.25/3 

WB Case B 0.003/0.25/3 

WC Case C 0.003/0.25/3 

SA Case A 0.1/90/920 

SB CaseB 0.1/90/920 

SC CaseC 0.1/90/920 



Note. — The left column lists the names of the different 2D simulations, middle column 
lists the method used to calculate the collision frequency between ions and neutrals. The 
last column shows the minimum, mean and maximum value of the unsigned magnetic field 
strength in the photosphere. 
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Fig. 18. — Following the Equation [161 the ratio between Pe^en and Pi^in is shown for the 
simulations labeled WB (top panel) and SB (bottom panel) at t = 500 s. The colorbar is 
the same for both panels and it is in logarithmic scale. 
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Fig. 19. — A study of the validity of the assumptions underlying the generalized Ohm's law. 
The dynamic frequency of the simulations 0.5 Hz) should remain lower than the frequency 
limits shown in the different panels for simulations WB (left panels) and SB (right panels) at 
t = 500 s. The frequencies are following the expressions Equation [T8l (top panels), Equation [3 
(middle panels) and Equation [T9l (bottom panels). The colorbar for each frequency is located 
at the top side and is in logarithmic scales. The white color is where the temperature is 
above 3 lO^K. 
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Table 2. Atomic info 



name 


H 


He 


C 


N 





Ne 


Na 


Mg 


abund 


12. 


11. 


8.55 


7.93 


8.77 


8.51 


6.18 


7.48 


mass ion 


1.008 


4.003 


12.01 


14.01 


16. 


20.18 


23. 


24.32 


Xi 


13.595 


24.58 


11.256 


14.529 


13.614 


21.559 


5.138 


7.644 


name 


Al 


Si 


S 


K 


Ca 


Cr 


Fe 


Ni 


abund 


6.4 


7.55 


7.21 


5.05 


6.33 


5.47 


7.5 


5.08 


mass ion 


26.97 


28.06 


32.06 


39.1 


40.08 


52.01 


55.85 


58.69 


Xi 


5.984 


8.149 


10.357 


4.339 


6.111 


6.763 


7.896 


7.633 



Note. — The atomic species, abundances (log of number of atoms per 10^^ Hydrogen 
atoms), mass ion (uma), and ionization fraction (eV) of the 16 most important atomic 
species in the solar atmosphere are listed from the top to the bottom row. The various 
collision frequencies and the electron density are calculated taking into account the atomic 
species in this table. 



